2.1 Histórico DHI

Para la elaboracion del modelo predictivo se cuenta con un Histórico de datos 15 minutales aportados por la presa de Belesar en la que se aportan datos como Aportación, Nivel del embalse, Caudal turbinado, lluvia y algunas otras variables.

DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_DHI_Belesar.RDS'))

DHI

plot(DHI$`APORTACION (m3/s)`,
     xlab = "Date",
     ylab= "Aportacion [m³/s]")

outlierKD <- function(dt, var) {
  var_name <- eval(substitute(var),eval(dt))
  tot <- sum(!is.na(var_name))
  na1 <- sum(is.na(var_name))
  m1 <- mean(var_name, na.rm = T)
  par(mar = rep(2, 4))
  par(mfrow=c(2, 2), oma=c(0,0,3,0))
  boxplot(var_name, main="With outliers")
  hist(var_name, main="With outliers", xlab=NA, ylab=NA)
  outlier <- boxplot.stats(var_name)$out
  mo <- mean(outlier)
  var_name <- ifelse(var_name %in% outlier, NA, var_name)
  boxplot(var_name, main="Without outliers")
  hist(var_name, main="Without outliers", xlab=NA, ylab=NA)
  title("Outlier Check", outer=TRUE)
  na2 <- sum(is.na(var_name))
  message("Outliers identified: ", na2 - na1, " from ", tot, " observations")
  message("Proportion (%) of outliers: ", (na2 - na1) / tot*100)
  message("Mean of the outliers: ", mo)
  m2 <- mean(var_name, na.rm = T)
  message("Mean without removing outliers: ", m1)
  message("Mean if we remove outliers: ", m2)
  
    dt[as.character(substitute(var))] <- invisible(var_name)
    assign(as.character(as.list(match.call())$dt), dt, envir = .GlobalEnv)
    message("Outliers successfully removed", "\n")
    return(invisible(dt))
   
}

DHI$`APORTACION (m3/s)`[which(DHI$`APORTACION (m3/s)`<0)]<- NA 

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 826 from 34552 observations

## Proportion (%) of outliers: 2.39059967585089

## Mean of the outliers: 1075.81016949153

## Mean without removing outliers: 131.482568013429

## Mean if we remove outliers: 108.354577773824

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 317 from 33726 observations

## Proportion (%) of outliers: 0.939927652256419

## Mean of the outliers: 410.450378548896

## Mean without removing outliers: 108.354577773824

## Mean if we remove outliers: 105.488153491574

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 228 from 33409 observations

## Proportion (%) of outliers: 0.682450836600916

## Mean of the outliers: 391.157280701754

## Mean without removing outliers: 105.488153491574

## Mean if we remove outliers: 103.525205991381

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 159 from 33181 observations

## Proportion (%) of outliers: 0.479189897833097

## Mean of the outliers: 380.017798742138

## Mean without removing outliers: 103.525205991381

## Mean if we remove outliers: 102.193901944158

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 168 from 33022 observations

## Proportion (%) of outliers: 0.5087517412634

## Mean of the outliers: 373.650892857143

## Mean without removing outliers: 102.193901944158

## Mean if we remove outliers: 100.805797771961

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 50 from 32854 observations

## Proportion (%) of outliers: 0.152188470201498

## Mean of the outliers: 365.8664

## Mean without removing outliers: 100.805797771961

## Mean if we remove outliers: 100.40179124497

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 4 from 32804 observations

## Proportion (%) of outliers: 0.0121936349225704

## Mean of the outliers: 362.2675

## Mean without removing outliers: 100.40179124497

## Mean if we remove outliers: 100.369856402439

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 1 from 32800 observations

## Proportion (%) of outliers: 0.00304878048780488

## Mean of the outliers: 361.84

## Mean without removing outliers: 100.369856402439

## Mean if we remove outliers: 100.361884508674

## Outliers successfully removed

outlierKD(DHI,`APORTACION (m3/s)`)

## Outliers identified: 0 from 32799 observations

## Proportion (%) of outliers: 0

## Mean of the outliers: NaN

## Mean without removing outliers: 100.361884508674

## Mean if we remove outliers: 100.361884508674

## Outliers successfully removed

plot(DHI$`APORTACION (m3/s)`,
     xlab = "Date",
     ylab= "Aportacion [m³/s]",
     type = "l")

sum(is.na(DHI$`APORTACION (m3/s)`))

## [1] 7180

library(imputeTS)

#Rellenamos Huecos
DHI$`APORTACION (m3/s)`<-  na.interpolation(DHI$`APORTACION (m3/s)`)


#Comprobamos
sum(is.na(DHI$`APORTACION (m3/s)`))

## [1] 0

library(TTR)

for(n in c(1,seq(48,192,length.out = 4))){
  Mavg<- SMA(DHI$`APORTACION (m3/s)`, n) 
  
  plot(DHI$DATE, Mavg,
       xlab = "Date",
       ylab= "Aportacion [m³/s]", 
       type = "l",
       main = paste0("Aportación en Belesar. MA= ", n))
}

outlierKD(DHI,`NIVEL EMBALSE (msnm)`)

## Outliers identified: 1 from 39978 observations

## Proportion (%) of outliers: 0.00250137575666617

## Mean of the outliers: 0

## Mean without removing outliers: 310.775323928161

## Mean if we remove outliers: 310.783097781224

## Outliers successfully removed

DHI$`NIVEL EMBALSE (msnm)`<-  na.interpolation(DHI$`NIVEL EMBALSE (msnm)`)

x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(c(0,diff(DHI$`NIVEL EMBALSE (msnm)`)*8000+150), x)

aportacion_SMA<- aportacion_SMA[!is.na(aportacion_SMA)]
nivel_SMA<- nivel_SMA[!is.na(nivel_SMA)]


plot(aportacion_SMA,type = "l",
     xlab =" " ,ylab = "",xaxt= 'n',yaxt='n',
     main=paste0("Variacion del nivel y aportacion\n Correlación de: ", 
                 round(cor(aportacion_SMA, nivel_SMA), 3)))
lines(nivel_SMA, col = "red")

ccf_belesar<- ccf(aportacion_SMA, nivel_SMA, lag.max = 5000)

#Máxima correlación
max(ccf_belesar$acf)

## [1] 0.5892262

#Con cuanto desfase se produce la máxima correlación. 
ccf_belesar$lag[which.max(ccf_belesar$acf)]

## [1] 34

#Retrasamos hacia detrás todos los datos de lluvia porque la acumulada de toda la hora la suma en la hora siguiente....
DHI$`LLUVIA ACUMULADA DÍA (l/m2)`<- lead(DHI$`LLUVIA ACUMULADA DÍA (l/m2)`)

Desacumular_lluvia_2018<- DHI[which(year(DHI$DATE)==2018),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
                                                                         lluvia=ifelse(desacumulada>=0, desacumulada, 0))


Desacumular_lluvia_2019<- DHI[which(year(DHI$DATE)==2019),] %>% group_by(yday(DATE)) %>% mutate(desacumulada= c(diff(`LLUVIA ACUMULADA DÍA (l/m2)`),0),
                                                                              lluvia=ifelse(desacumulada>=0, desacumulada, 0))
Desacumular_lluvia<- rbind(Desacumular_lluvia_2018, Desacumular_lluvia_2019)


Medias_horarias_2018<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2018,]) %>% group_by(yday(DATE), hour(DATE))  %>% 
  summarize(., Acum_horaria=sum(lluvia, na.rm = T),
            aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
            nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))

Medias_horarias_2019<- as.data.frame(Desacumular_lluvia[year(Desacumular_lluvia$DATE)==2019,]) %>% group_by(yday(DATE),hour(DATE))  %>% 
  summarize(., Acum_horaria=sum(lluvia, na.rm = T),
            aport_mean=mean(`APORTACION (m3/s)`, na.rm = T),
            nivel_mean=mean(`NIVEL EMBALSE (msnm)`, na.rm = T))

Medias_horarias<- rbind(Medias_horarias_2018, Medias_horarias_2019)
vector_Date<- seq(range(DHI$DATE)[1],range(DHI$DATE)[2],
                  by="hour")

Lluvia_acum_horaria<- as.data.frame(cbind(as.character(vector_Date[2:length(vector_Date)]), 
                                          Medias_horarias[,3:5]))



Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`<- ymd_hms(Lluvia_acum_horaria$`as.character(vector_Date[2:length(vector_Date)])`)


colnames(Lluvia_acum_horaria)<- c("Date", "Lluvia_mm", "aport_mean", "nivel_mean")

for (n in seq(6,24*3, by=6)) {
  y<- n
  aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
  nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)
  
  
  
  Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria, 
                                  
                                  aportacion_mean_SMA,
                                  nivel_mean_SMA))
  
  
  
  colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_mean_SMA",  "nivel_mean_SMA")
  Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]
  
  
  
  ccf_aport_mean<- ccf(Tabla_DHI$aport_mean_SMA,
                       Tabla_DHI$nivel_mean_SMA, 
                       lag.max = 5000, 
                       plot = T)
  
  print(paste0("Máxima correlacion de ",
               round(max(ccf_aport_mean$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_aport_mean$lag[which.max(ccf_aport_mean$acf)]/24, digits = 2), 
               " días. Con un SMA de :", 
               n, 
               " horas"))
  
}

## [1] "Máxima correlacion de 0.609 .Para un desfase de: -76.79 días. Con un SMA de :6 horas"

## [1] "Máxima correlacion de 0.616 .Para un desfase de: -76.58 días. Con un SMA de :12 horas"

## [1] "Máxima correlacion de 0.621 .Para un desfase de: -76.42 días. Con un SMA de :18 horas"

## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días. Con un SMA de :24 horas"

## [1] "Máxima correlacion de 0.627 .Para un desfase de: -76.29 días. Con un SMA de :30 horas"

## [1] "Máxima correlacion de 0.629 .Para un desfase de: -76.25 días. Con un SMA de :36 horas"

## [1] "Máxima correlacion de 0.632 .Para un desfase de: -76.21 días. Con un SMA de :42 horas"

## [1] "Máxima correlacion de 0.634 .Para un desfase de: -76.08 días. Con un SMA de :48 horas"

## [1] "Máxima correlacion de 0.636 .Para un desfase de: -75.92 días. Con un SMA de :54 horas"

## [1] "Máxima correlacion de 0.638 .Para un desfase de: -75.88 días. Con un SMA de :60 horas"

## [1] "Máxima correlacion de 0.64 .Para un desfase de: -75.79 días. Con un SMA de :66 horas"

## [1] "Máxima correlacion de 0.642 .Para un desfase de: -75.67 días. Con un SMA de :72 horas"

#se selecciona 96 porque es la cantidad de datos que corresponden a 1 día. 
x<- 96
aportacion_SMA<- SMA(DHI$`APORTACION (m3/s)`,x)
nivel_SMA<- SMA(DHI$`NIVEL EMBALSE (msnm)`, x)

aportacion_horaria<- aportacion_SMA[DHI$DATE%in%vector_Date]
nivel_horario<- nivel_SMA[DHI$DATE%in%vector_Date]

#Se ponen 24 horas para que equivalgan a 1 día 
y<- 24
aportacion_mean_SMA<- SMA(Lluvia_acum_horaria$aport_mean,y)
nivel_mean_SMA<- SMA(Lluvia_acum_horaria$nivel_mean,y)





Tabla_DHI<- as.data.frame(cbind(Lluvia_acum_horaria, 
                                aportacion_horaria[2:length(aportacion_horaria)],
                                nivel_horario[2:length(nivel_horario)], 
                                aportacion_mean_SMA,
                                nivel_mean_SMA))



colnames(Tabla_DHI)<- c(names(Lluvia_acum_horaria), "aport_SMA", "nivel_SMA", 
                        "aport_mean_SMA",  "nivel_mean_SMA")

Tabla_DHI<- Tabla_DHI[complete.cases(Tabla_DHI),]

 k<- max(Tabla_DHI$nivel_mean_SMA)/max(Tabla_DHI$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_mean_SMA/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_mean_SMA),
                                                        max(Tabla_DHI$nivel_mean_SMA),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean_SMA),
                                    max(Tabla_DHI$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Nivel y aportación. SMA (24 horas) de la media horaria")

ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA, 
    Tabla_DHI$nivel_mean_SMA, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.583 .Para un desfase de: -41.67 días"

 k<- max(Tabla_DHI$nivel_SMA)/max(Tabla_DHI$aport_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_SMA/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_SMA),
                                                        max(Tabla_DHI$nivel_SMA),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_SMA),
                                    max(Tabla_DHI$aport_SMA),
                                    by=20)) + 
    ggtitle("Nivel y aportación. SMA (96 datos = 1 día) de los datos 15-minutales")

ccf_DHI<- ccf(Tabla_DHI$aport_SMA, 
    Tabla_DHI$nivel_SMA, 
    lag.max = 4000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.624 .Para un desfase de: -76.38 días"

 k<- max(Tabla_DHI$nivel_mean)/max(Tabla_DHI$aport_mean)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=nivel_mean/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI$nivel_mean),
                                                        max(Tabla_DHI$nivel_mean),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean),
                                    max(Tabla_DHI$aport_mean),
                                    by=20)) + 
    ggtitle("Nivel y aportación. Medias horarias")

ccf_DHI<- ccf(Tabla_DHI$aport_mean, 
    Tabla_DHI$nivel_mean, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.54 .Para un desfase de: -41.67 días"

 k<- max(Tabla_DHI$Lluvia_mm)/max(Tabla_DHI$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Lluvia Horaria [l/m²]", 
                                           breaks = seq(min(Tabla_DHI$Lluvia_mm),
                                                        max(Tabla_DHI$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI$aport_mean_SMA),
                                    max(Tabla_DHI$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion")

ccf_DHI<- ccf(Tabla_DHI$aport_mean_SMA, 
    Tabla_DHI$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.109 .Para un desfase de: 1.71 días"

Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12), ]
 k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI_cut, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
                                                        max(Tabla_DHI_cut$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
                                    max(Tabla_DHI_cut$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion. Diciembre")

ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA, 
    Tabla_DHI_cut$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.371 .Para un desfase de: 1.62 días"

Tabla_DHI_cut<- Tabla_DHI[which(month(Tabla_DHI$Date)==12 | year(Tabla_DHI$Date)==2019), ]
 k<- max(Tabla_DHI_cut$Lluvia_mm)/max(Tabla_DHI_cut$aport_mean_SMA)
 
  ggplot(data=Tabla_DHI_cut, aes(x=Date))+
    geom_line(aes(y=aport_mean_SMA))+
    xlab("Date")+ylab("Aportación [m³/s]")+
    theme(panel.background = element_blank(),
          panel.grid = element_blank())+
    geom_line(aes(y=Lluvia_mm/k),group=1, color="red") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k), name = "Nivel [msnm]", 
                                           breaks = seq(min(Tabla_DHI_cut$Lluvia_mm),
                                                        max(Tabla_DHI_cut$Lluvia_mm),
                                                        by=10)),
                       breaks = seq(min(Tabla_DHI_cut$aport_mean_SMA),
                                    max(Tabla_DHI_cut$aport_mean_SMA),
                                    by=20)) + 
    ggtitle("Lluvia y aportacion. Diciembre-Enero-Febrero")

ccf_DHI<- ccf(Tabla_DHI_cut$aport_mean_SMA, 
    Tabla_DHI_cut$Lluvia_mm, 
    lag.max = 1000)

 print(paste0("Máxima correlacion de ",
               round(max(ccf_DHI$acf), digits = 3) ,
               " .Para un desfase de: ",
               round(ccf_DHI$lag[which.max(ccf_DHI$acf)]/24, digits = 2), 
        " días"))

## [1] "Máxima correlacion de 0.332 .Para un desfase de: 1.71 días"

path_dataDHI <- here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS')

#saveRDS(Tabla_DHI,path_dataDHI)

2.2 Datos WRF

Los datos del histórico de WRF están organizados en formato lista… hay para cada localización 2 dataframes;

• D1: WRF predicción del mismo día
• D2: WRF predicción de mañana

A continuación se presenta el loop y las funciones que se emplean para generar todo el histórico WRF de Belesar… Cabe mencionar que contamos con una carpeta con Archivos RDS ordenados por fechas para la localización de belesar.

Usamos las siguientes funciones.

lon_lat_df_ls<- function(parque_list){
  fecha<- names(parque_list)
  fecha[!str_detect(fecha, " ")]<- paste0(fecha[!str_detect(fecha, " ")], " 00:00:00")
  
  n_lon<- parque_list[[1]]$lon
  n_lat<- parque_list[[1]]$lat
  
  nombres<- paste0(n_lon,"__",n_lat)
  
  
  list_localizaciones<- list()
  for (localizaciones in 1:length(parque_list[[1]][,1]) ) {
    l<- lapply(parque_list, function(x) x[localizaciones,])
    df <- data.frame(matrix(unlist(l), nrow=length(l), byrow=T))
    df<- as.data.frame(cbind(ymd_hms(fecha), df))
    colnames(df)<-c("fechas", colnames(parque_list[[1]]))
    list_localizaciones[[localizaciones]]<- df
  }
  
  names(list_localizaciones)<- nombres
  return(list_localizaciones)
  
  
}
uv_transformation<- function(tabla_comp){
  
  nombres<- colnames(tabla_comp)
  
  u10<- tabla_comp$U10_MEAN
  v10<- tabla_comp$V10_MEAN
  
  u10_max<- tabla_comp$U10_MAX
  v10_max<- tabla_comp$V10_MAX
  
  wind_abs <- sqrt(u10^2 + v10^2)
  wind_dir_rad <-  atan2(u10/wind_abs, v10/wind_abs) 
  wind_dir_deg1 <-  wind_dir_rad * 180/pi 
  wind_dir_deg2 <-  wind_dir_deg1+ 180 
  
  
  wind_abs_max <- sqrt(u10_max^2 + v10_max^2)
  wind_dir_rad_max <-  atan2(u10_max/wind_abs_max, v10_max/wind_abs_max) 
  wind_dir_deg1_max <-  wind_dir_rad_max * 180/pi 
  wind_dir_deg2_max <-  wind_dir_deg1_max+ 180 
  
  tabla_comp<- as.data.frame(cbind(tabla_comp,wind_abs,wind_dir_deg2,
                                   wind_abs_max,wind_dir_deg2_max))
  colnames(tabla_comp)<- c(nombres,"WS","WD","WS_MAX","WD_MAX")
  tabla_comp$U10_MAX<- NULL
  tabla_comp$V10_MAX<- NULL
  tabla_comp$U10_MEAN<- NULL
  tabla_comp$V10_MEAN<- NULL
  return(tabla_comp)
  
}
extract_rain_data<- function(Belesar_lolat_df){
  
  rain<- Belesar_lolat_df[,c(1,2,3,4,5)]
  rain$pre_acum<- rain$RAINC+rain$RAINNC
  rain$RAINC<- NULL
  rain$RAINNC<- NULL
  
  prep_hourly<- vector()
  for (i in 1:length(rain$pre_acum)) {
    if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
      prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
    }
    
  }
  rain$prep_hourly<- prep_hourly 
  
  return(rain)
}
extract_rain_data2<- function(Belesar_lolat_df){
  
  rain<- Belesar_lolat_df[,c(1,2,3,4,5,6,10,17,20)]
  rain$pre_acum<- rain$RAINC+rain$RAINNC
  
  prep_hourly<- vector()
  for (i in 1:length(rain$pre_acum)) {
    if(i==1){prep_hourly[i]<- rain$pre_acum[i]}else{
      prep_hourly[i]<- rain$pre_acum[i]-rain$pre_acum[i-1]
    }
    
  }
  rain$prep_hourly<- prep_hourly 
  
  return(rain)
}
Return_periodo_Belesar<- function(){
  All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'), 
                                 recursive = F, full.names = T)
  RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]
  
  RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]
  
  
  Belesar_data<- readRDS(RDS_Belesar1[1])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
  
  fecha_ini<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[1]
  
  Belesar_data<- readRDS(RDS_Belesar1[length(RDS_Belesar)])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data)
  
  fecha_last<- Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas[length(Belesar_rain$`-8.02328491210938__42.1343421936035`$fechas)]
  
  periodo_WRF<- seq(fecha_ini, fecha_last, by="hour")
  
  Tabla_WRF<- as.data.frame(matrix(ncol = 10, nrow = length(periodo_WRF)))
  colnames(Tabla_WRF)<- colnames(Belesar_rain[[1]])
  Tabla_WRF$fechas<- periodo_WRF
  return(Tabla_WRF)
}

#Listamos archivos dentro de la carpeta de Belesar
All_files_Belesar<- list.files(here::here('Data/Parques/Belesar/'),
                               recursive = F, full.names = T)


#Detectamos cuales son RDS
RDS_Belesar<- All_files_Belesar[str_detect(All_files_Belesar, ".RDS")]

#Eliminamos los RDS que no son de WRF
RDS_Belesar1<- RDS_Belesar[!str_detect(RDS_Belesar, "E001")]

#Construimos Lista para cada instante de tiempo
Lista_localizacion<- list() 
for (i in 1:length(RDS_Belesar1)) {
  Belesar_data<- readRDS(RDS_Belesar1[i])
  Belesar_lolat<- lon_lat_df_ls(Belesar_data)
  Belesar_lolat1<- lapply(Belesar_lolat, uv_transformation)
  Belesar_rain<- lapply(Belesar_lolat1, extract_rain_data2)
  Lista_localizacion[[i]]<- Belesar_rain
}

#Nombramos la lista
names_fechas<- sapply(RDS_Belesar1, function(x){
  r<- str_split(x, "/")
  return(str_remove(str_remove(r[[1]][length(r[[1]])], ".RDS"), "Belesar_"))
})
names(Lista_localizacion)<- names_fechas


#Creamos una lista por localización
Lista_localizacion2<- list()
for (i in 1:length(Lista_localizacion[[1]])) {
  Lista_localizacion2[[i]]<- lapply(Lista_localizacion, 
                               function(x) return(x[[i]]))
}

#Nombramos la lista
names(Lista_localizacion2)<- names(Lista_localizacion[[1]])

#Guardamos la lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
#saveRDS(Lista_localizacion2, nombre_lista)

#Cargamos lista
path_lista_total<- here::here('Data/Parques/Belesar/Historico/')
nombre_lista<- paste0(path_lista_total, 'Historico_WRF_Belesar_Variables.RDS')
Lista_total1<- readRDS(nombre_lista)

#Juntamos todos los Dataframes
Lista_total_MF<- lapply(Lista_total1, function(x) bind_rows(x))

#creamos dos data.frames... uno para D1 y otro para D2
Lista_d1_d2_loc<- list()
for (i in 1:length(Lista_total_MF)) {
  p<- Lista_total_MF[[i]]
  d1<- p[duplicated(p$fechas),]
  d1<-d1[!duplicated(d1$fechas),]
  d2<- p[!duplicated(p$fechas),]
  
  d2_qneed1<-d2[!(d2$fechas%in%d1$fechas),]
  
  
  d1_2<-bind_rows(d1,d2_qneed1)
  d1_2<-d1_2[order(d1_2$fechas),]
  
  d2<-d2[order(d2$fechas),]
  
  d1_2$pre_acum<- NULL
  d2$pre_acum<- NULL
  
  colnames(d1_2)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
  
  colnames(d2)<-  c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")
  
  lista_loc_d12<- list(d1_2,d2)
  names(lista_loc_d12)<- c("D1", "D2")
  
  Lista_d1_d2_loc[[i]]<- lista_loc_d12
}

#nombramos la lsta
names(Lista_d1_d2_loc)<- names(Lista_total_MF)

Tabla_periodo1<- Return_periodo_Belesar()
colnames(Tabla_periodo1)<- c("Date", "LON", "LAT", "RAINC", "RAINNC","RAINSH", "T02_MEAN","PSFC","WS_MAX", "prep_hourly")

Lista_d1_d2_loc2<- list()
for (j in 1:length(Lista_d1_d2_loc)) {
  prueba_list<- Lista_d1_d2_loc[[j]]
  lista_retorno<- list()
  for(i in 1:2){
    prueba<- prueba_list[[i]]
    Tabla_periodo<- Tabla_periodo1
    Tabla_periodo$LON[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LON
    Tabla_periodo$LAT[match(prueba$Date,Tabla_periodo$Date)] <- prueba$LAT
    Tabla_periodo$RAINC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINC
    Tabla_periodo$RAINNC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINNC
    Tabla_periodo$RAINSH[match(prueba$Date,Tabla_periodo$Date)] <- prueba$RAINSH
    Tabla_periodo$T02_MEAN[match(prueba$Date,Tabla_periodo$Date)] <- prueba$T02_MEAN
    Tabla_periodo$PSFC[match(prueba$Date,Tabla_periodo$Date)] <- prueba$PSFC
    Tabla_periodo$WS_MAX[match(prueba$Date,Tabla_periodo$Date)] <- prueba$WS_MAX
    Tabla_periodo$prep_hourly[match(prueba$Date,Tabla_periodo$Date)] <- prueba$prep_hourly
    lista_retorno[[i]]<- Tabla_periodo
  }
  names(lista_retorno)<- c("D1", "D2")
  Lista_d1_d2_loc2[[j]]<- lista_retorno 
  
  
}
names(Lista_d1_d2_loc2)<- names(Lista_d1_d2_loc)

#Creamos lista con las variables afinadas
Lista_d1_d2_loc3<- lapply(Lista_d1_d2_loc2, function(x){
  x[[1]]$RAINC<- NULL
  x[[1]]$RAINNC<- NULL
  x[[1]]$RAINSH<- NULL
  x[[2]]$RAINC<- NULL
  x[[2]]$RAINNC<- NULL
  x[[2]]$RAINSH<- NULL
  
  return(x)})
  
names(Lista_d1_d2_loc3)<- names(Lista_d1_d2_loc2)

#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS')
#saveRDS(Lista_d1_d2_loc3,path_hist_WRF)

2.3 Juntando históricos

Juntamos a cada data.frame de WRF la tabla afinada de DHI.

#Usamos dplyr para juntar ambos datasets

Belesar_WRF<- readRDS(here::here('Data/Parques/Belesar/Historico/Historico_WRF_Belesar_Variables_D1D2.RDS'))
Belesar_DHI<- readRDS(here::here('Data/Parques/Belesar/Historico/DHI_historico_afinado.RDS'))

df2<- Belesar_DHI

Belesar_Merge<- list()
for (j in 1:length(Belesar_WRF)) {
  lista_retorno<- list()
  for(i in 1:2){
    df1<-  Belesar_WRF[[j]][[i]]
    Merge_table<- left_join(df1, df2, by=c("Date"))
    lista_retorno[[i]]<- Merge_table
  }
  names(lista_retorno)<- c("D1", "D2")
  Belesar_Merge[[j]]<- lista_retorno 
}
names(Belesar_Merge)<- names(Belesar_WRF)

#Belesar merge completecases
Belesar_Merge_cc<- list()
for (j in 1:length(Belesar_Merge)) {
    lista_retorno<- list()
  for(i in 1:2){
    df1<-  Belesar_Merge[[j]][[i]]
    df1$Acumulated<- NULL
    Table_fine<- df1[complete.cases(df1),]
    lista_retorno[[i]]<- Table_fine
  }
  names(lista_retorno)<- c("D1", "D2")
  Belesar_Merge_cc[[j]]<- lista_retorno 
}
names(Belesar_Merge_cc)<- names(Belesar_Merge)



#Guardamos
path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
#saveRDS(Belesar_Merge_cc,path_hist_WRF)

4.1 Octubre-Febrero

#Empleamos los datos a partir de octubre y eliminamos datos negativos de 
# precipitacion WRF que hemos visto que hay alguno. 
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) {
    r<- r[which(r$Date > ymd("2018/10/01")),]
    r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly) 
    return(r)
  })
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_oct)) {
  Corr<- cor(clean_data_oct[[i]][[1]]$prep_hourly, 
      clean_data_oct[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_oct)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
                      c(as.numeric(LON),
                        as.numeric(LAT)),
                      fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}

Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order

#Plot_rain3 necesita columna Date y columna rain
plot_rain3<- function(data,data2, titulo){
  ggplot(data=data, aes(x=Date))+
    geom_line(aes(y=rain), stat="identity")+
    xlab("Date")+ylab("Lluvia por hora [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank()) +
    geom_line(data= data2, aes(x=Date, y=rain), color="red", alpha=0.5)+
    geom_text(aes(as.POSIXct(ymd("2019/02/01")),10, label=paste0("Cor: ",round(cor(data$rain,data2$rain), 
                                                                               digits = 2)))) +
    ggtitle(titulo) + 
     theme(plot.title = element_text(hjust = 0.5))
  
}

#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
                              clean_data_oct[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_oct[[i]][[1]]$Date),
                              clean_data_oct[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_oct)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

4.2 Enero-Febrero

Es reseñable que disminuyendo el tiempo de ploteo la correlación aumenta.

clean_data_jan<-lapply(clean_data_oct, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/01")),])
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
  Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly, 
             clean_data_jan[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}

Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]

Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]

Cor_order

#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

4.3 15 de Enero- 10 de febrero

Apuramos más el rango para que únicamente coja un periodo lluvioso como fue el de el 15 de enero al pocos de febrero

clean_data_jan<-lapply(clean_data_oct, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2019/01/15") & r$Date < ymd("2019/02/10")),])
})

Cor_rain_place<- data.frame(matrix(ncol = 4))
colnames(Cor_rain_place)<- c("LON", "LAT", "Corr", "Dist")
for (i in 1:length(clean_data_jan)) {
  Corr<- cor(clean_data_jan[[i]][[1]]$prep_hourly, 
             clean_data_jan[[i]][[1]]$Lluvia_mm)
  LON<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(clean_data_jan)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, Corr, Dist))
  
  
}

Cor_order<- Cor_rain_place[order(Cor_rain_place$Corr, decreasing = T),]
Cor_order

#Ploteamos los 5 puntos con mejor correlacion
for (i in order(Cor_rain_place$Corr, decreasing = T)[1:5]) {
  data1<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$prep_hourly))
  
  names(data1)<- c("Date", "rain")
  data1$Date<- ymd_hms(data1$Date)
  data1$rain<- as.numeric(as.character(data1$rain))
  
  
  
  data2<- as.data.frame(cbind(as.character(clean_data_jan[[i]][[1]]$Date),
                              clean_data_jan[[i]][[1]]$Lluvia_mm))
  
  names(data2)<- c("Date", "rain")
  data2$Date<- ymd_hms(data2$Date)
  data2$rain<- as.numeric(as.character(data2$rain))
  
  
  x<- plot_rain3(data = data1 ,
             data2 = data2,
             titulo = paste0(paste(round(as.numeric(str_split(names(clean_data_jan)[i],
                                                       "_")[[1]][c(1,3)]),
                                  digits = 3), 
                            collapse = " "), "\n n = ", 
                            i, " Dist = ",
                            round((Cor_rain_place$Dist[i]/1000), digits = 2), " km" ))
  print(x)
  
  
  
  
}

5.1 Reducimos el tiempo

Antes estabamos usando los datos desde octubre hasta febrero… Ahora empleamos solo los datos desde diciembre hasta febrero.

clean_data_dec<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2018/12/01")),])
})
Daiyli_avgs<- lapply(clean_data_dec, function(x){
  x[[1]] %>% group_by(year(Date), yday(Date)) %>%  
    summarize(rWRF_acum= sum(prep_hourly),
              rDHI_acum= sum(Lluvia_mm),
              rWRF_mean= sum(prep_hourly),
              rDHI_mean= sum(Lluvia_mm),
              aport_SMA_diario=mean(aport_SMA),
              aport_mean_diario=mean(aport_mean),
              aport_mean_SMA_diario=mean(aport_mean_SMA),
              nivel_SMA_diario=mean(nivel_SMA),
              nivel_mean_diario=mean(nivel_mean),
              nivel_mean_SMA_diario=mean(nivel_mean_SMA))

  
  })

Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
  
  x$Date<- ymd(ifelse(x$`year(Date)`==2018,
                      as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
                      as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
  x$`year(Date)`<- NULL
  x$`yday(Date)`<- NULL
  return(x)
})

Cor_rain_place<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  Corr1<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
             Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr2<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr3<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr4<- cor(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum)
  
  Corr5<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario)
  Corr6<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario)
  Corr7<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario)
  Corr8<- cor(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place[i,]<- as.data.frame(cbind(LON, LAT, round(Corr1, digits = 2),
                                           round(Corr2, digits = 2),
                                           round(Corr3, digits = 2),
                                           round(Corr4, digits = 2),
                                           round(Corr5, digits = 2),
                                           round(Corr6, digits = 2),
                                           round(Corr7, digits = 2),
                                           round(Corr, digits = 2),
                                           Dist))
  
  
}

colnames(Cor_rain_place)<- c("LON","LAT","AWRF_aport1",
                             "AWRF_aport2", "AWRF_aport3",
                             "AWRF_ADHI","MWRF_aport1",
                             "MWRF_aport2", "MWRF_aport3",
                             "MWRF_MDHI","Dist")

Cor_order<- Cor_rain_place[order(Cor_rain_place$Dist),]
Cor_order

Cor_rain_place_ccf<- data.frame(matrix(ncol = 11))
for (i in 1:length(Daiyli_avgs2)) {
  ccfr1<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr2<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr3<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr4<- ccf(Daiyli_avgs2[[i]]$rWRF_acum, 
              Daiyli_avgs2[[i]]$rDHI_acum, plot = FALSE)
  
  ccfr5<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_SMA_diario, plot = FALSE)
  ccfr6<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_diario, plot = FALSE)
  ccfr7<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$aport_mean_SMA_diario, plot = FALSE)
  ccfr8<- ccf(Daiyli_avgs2[[i]]$rWRF_mean, 
              Daiyli_avgs2[[i]]$rDHI_mean, plot = FALSE)
  
  LON<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][1])
  LAT<- as.numeric(str_split(names(Daiyli_avgs2)[i], "_")[[1]][3])
  Dist<- distm(c(-7.713948, 42.628577), 
               c(as.numeric(LON),
                 as.numeric(LAT)),
               fun = distHaversine)
  Cor_rain_place_ccf[i,]<- as.data.frame(cbind(LON, LAT, 
                                               round(max(ccfr1$acf), digits = 2),
                                               round(max(ccfr2$acf), digits = 2),
                                               round(max(ccfr3$acf), digits = 2),
                                               round(max(ccfr4$acf), digits = 2),
                                               round(max(ccfr5$acf), digits = 2),
                                               round(max(ccfr6$acf), digits = 2),
                                               round(max(ccfr7$acf), digits = 2),
                                               round(max(ccfr8$acf), digits = 2), Dist))
  print(paste("Mejor corr para un retardo de: ", 
              ccfr1$lag[which.max(ccfr1$acf)],
              ccfr2$lag[which.max(ccfr2$acf)],
              ccfr3$lag[which.max(ccfr3$acf)],
              ccfr4$lag[which.max(ccfr4$acf)],
              ccfr5$lag[which.max(ccfr5$acf)],
              ccfr6$lag[which.max(ccfr6$acf)],
              ccfr7$lag[which.max(ccfr7$acf)],
              ccfr8$lag[which.max(ccfr8$acf)], sep = " "))  
  
}

## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -6 -5 -6 0 -6 -5 -6 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -4 -3 0 -3 -4 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -2 -2 -2 0 -2 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -2 0 -3 -2 -2 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -3 -2 -3 0 -3 -2 -3 0"
## [1] "Mejor corr para un retardo de:  -13 -12 -13 0 -13 -12 -13 0"

colnames(Cor_rain_place_ccf)<-  c("LON","LAT","AWRF_aport1",
                                  "AWRF_aport2", "AWRF_aport3",
                                  "AWRF_ADHI","MWRF_aport1",
                                  "MWRF_aport2","MWRF_aport3",
                                  "MWRF_MDHI","Dist")

Cor_order<- Cor_rain_place_ccf[order(Cor_rain_place_ccf$Dist),]
Cor_order

A continuación extraemos la media y el máximo de cada columna.

#Media de cada columna
colMeans(Cor_rain_place_ccf[3:(length(Cor_rain_place_ccf)-1)])

## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##   0.2883721   0.2939535   0.2890698   0.6567442   0.2883721   0.2939535 
## MWRF_aport3   MWRF_MDHI 
##   0.2890698   0.6567442

#función ColMAX
colMax <- function(data) sapply(data, max, na.rm = TRUE)

colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]

## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##        0.49        0.50        0.49        0.73        0.49        0.50 
## MWRF_aport3   MWRF_MDHI 
##        0.49        0.73

De aquí podemos extraer las siguientes conclusiones:

• No hay diferencia en la correlación usando la lluvia acumulada o la lluvia media diaria.
• Para realizar correlación lluvia WRF y aportación el mejor metodo es el 3, es decir, haciendo SMA sobre los datos medios horarios de aportación. A partir de ahora emplearemos AWRF_aport3 y AWRF_ADHI, como indicadores para la correlacion lluvia WRF-aportación y lluviaWRF- lluvia observada
• Existe un punto con una correlación entre lluvia registrada y lluvia WRF de 92%. BUENAS NOTICIAS. Esto nos da información de que el modelo está haciendo bien su trabajo.

#función ColMAX
which_colMax <- function(data) sapply(data, which.max)

which_colMax(Cor_rain_place_ccf)[3:(length(Cor_rain_place_ccf)-1)]

## AWRF_aport1 AWRF_aport2 AWRF_aport3   AWRF_ADHI MWRF_aport1 MWRF_aport2 
##          18          18          18           1          18          18 
## MWRF_aport3   MWRF_MDHI 
##          18           1

Vemos que en cuanto a correlación entre WRF y aportación el mejor punto es el numero 40.

Cor_rain_place_ccf[40,]

Apreciamos que se encuentra a una distancia de 82 km. Esto es peligroso.

En cuanto a la lluvia…

Cor_rain_place_ccf[4,]

Vemos que se encuentra a 52 km.

Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T),c("AWRF_aport3","AWRF_ADHI", "Dist")]

Por lo visto en la tabla anterior, habrá que buscar un punto que tenga buenas correlaciones para ambas comparaciones y además esté relativamente cerca. En mi opinión estas son las mejores opciones. No sé cual elegir.

Cor_rain_place_ccf[c(19,23,24,14),c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]

Construimos una lista solo con los puntos las variables que nos interesan. Hemos ido reduciendo el rango de tiempo viendo como las correlaciones aumentaban, ahora que hemos escogido el mejor punto vamos a cojer los datos desde finales de octubre.

clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) {
    r<- r[which(r$Date > ymd("2018/10/01")),]
    r$prep_hourly<- ifelse(r$prep_hourly < 0,0,r$prep_hourly) 
    return(r)
  })
})
Daiyli_avgs<- lapply(clean_data_oct, function(x){
  x[[1]] %>% group_by(year(Date), yday(Date)) %>%  
    summarize(rWRF_acum= sum(prep_hourly),
              rDHI_acum= sum(Lluvia_mm),
              rWRF_mean= sum(prep_hourly),
              rDHI_mean= sum(Lluvia_mm),
              aport_SMA_diario=mean(aport_SMA),
              aport_mean_diario=mean(aport_mean),
              aport_mean_SMA_diario=mean(aport_mean_SMA),
              nivel_SMA_diario=mean(nivel_SMA),
              nivel_mean_diario=mean(nivel_mean),
              nivel_mean_SMA_diario=mean(nivel_mean_SMA))

  
  })

Daiyli_avgs2<- lapply(Daiyli_avgs, function(x){
  
  x$Date<- ymd(ifelse(x$`year(Date)`==2018,
                      as.character(as.Date(x$`yday(Date)`, origin= "2018-01-01")),
                      as.character(as.Date(x$`yday(Date)`, origin= "2019-01-01"))))
  x$`year(Date)`<- NULL
  x$`yday(Date)`<- NULL
  return(x)
})

Daiyli_avgs3<-lapply(Daiyli_avgs2[c(19,23,24,14)], function(x)x[, c("Date","rWRF_acum",
                                      "aport_mean_SMA_diario", 
                                      "rDHI_acum")])


head(Daiyli_avgs3[[1]])

range(Daiyli_avgs3[[1]]$Date)

## [1] "2018-10-26" "2019-02-22"

Apreciamos que tenemos un un dato diario entre el 26 de octubre y el 22 de febrero.

lapply(Daiyli_avgs3, function(x, titulo){
   k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
  ggplot(data=x, aes(x=Date))+
    geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
    geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+
    theme(panel.background = element_blank(),panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
                                           name = " LLuvia acumulada [mm]", 
                                           breaks = seq(min(cumsum(x$rDHI_acum)),
                                                        max(cumsum(x$rDHI_acum)),
                                                        by=50)),
                       breaks = seq(min(x$rDHI_acum),
                                    max(x$rDHI_acum),
                                    by=10))+
    ggtitle(titulo) +
    theme(plot.title = element_text(hjust = 0.5))

}, titulo=names(Daiyli_avgs3))

## $`-8.05770874023438__42.6729469299316`

## 
## $`-7.32501220703125__42.6962661743164`

## 
## $`-7.88290405273438__42.8138198852539`

## 
## $`-8.04904174804688__42.5383224487305`

Es interesante apreciar como pese a las altas correlaciones existe una pequeña diferencia al final del periodo en lo que a lluvia acumulada se refiere. Esto nos da pie a investigar cuales son los puntos que calculan mejor la lluvia acumulada.

vec_diff<- as.vector(unlist(lapply(Daiyli_avgs2, function(x){
  r<- sum(x$rWRF_acum)-sum(x$rDHI_acum)
  return(r)
})))

vec_diff[order(abs(vec_diff))]

##  [1]    4.404910    4.729176  -10.388001  -20.714521   31.905542
##  [6]  -32.582030   36.890532   39.813112  -49.308699  -49.663263
## [11]  -55.735768   76.225877  -77.398381  -79.743920  -81.969871
## [16]  -87.654145  -91.766186   91.816973  -98.120085  102.573777
## [21] -108.333345 -108.674158 -111.758937 -117.214062 -120.455026
## [26] -127.220242 -130.881619 -135.337830 -136.514876 -147.815890
## [31]  149.886050 -155.839446 -160.733827 -171.574046  180.299144
## [36] -185.876094  191.348111 -193.420657 -213.624562 -220.962254
## [41] -225.314246 -229.023156 -256.119342

Cor_rain_place_ccf[order(abs(vec_diff))[1:10],c("LON","LAT","AWRF_aport3","AWRF_ADHI", "Dist")]

Vemos que entre las dos tablas que hemos sacado… la localización que se repite es la número 14, que se encuentra a 29 km. Esto es muy interesante para realizar luego las predicciones en base a esta localización.

Cor_rain_place_ccf[order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10],c("AWRF_aport3","AWRF_ADHI", "Dist")]

Vemos que los puntos que se repiten en ambas listas… buscando mejor correlación con aportación, cercania y mejor resultado con la lluvia acumulada.

best_points<- order(abs(vec_diff))[1:10][order(abs(vec_diff))[1:10]%in%order(Cor_rain_place_ccf$AWRF_aport3, decreasing = T)[1:10]]

best_points

## [1] 28 18 19 32

litros_acum_diff<- vec_diff[best_points]
cbind(litros_acum_diff,
Cor_rain_place_ccf[best_points,c("AWRF_aport3","AWRF_ADHI", "Dist")])

Para un futuro tendremos que elegir entre estos 3 puntos. Añadimos también datos de correlación.

Daiyli_avgs4<-lapply(Daiyli_avgs2[c(14,18,19)], function(x)x[, c("Date","rWRF_acum",
                                      "aport_mean_SMA_diario", 
                                      "rDHI_acum")])

lapply(Daiyli_avgs4, function(x, titulo){
   k<- (max(c(cumsum(x$rWRF_acum),cumsum(x$rWRF_acum)))/max(c(x$rWRF_acum, x$rDHI_acum)))
  ggplot(data=x, aes(x=ymd(Date)))+
    geom_bar(aes(y=rWRF_acum), stat="identity", fill="red", alpha=0.5, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+theme(panel.background = element_blank(), 
                                                      panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rWRF_acum)/k), group = 1, col="red") +
    geom_bar(aes(y=rDHI_acum), stat="identity", fill="blue", alpha=0.4, size=0.05)+
    xlab("Date")+ylab("Lluvia diaria [mm/h]")+
    theme(panel.background = element_blank(),panel.grid = element_blank())+
    geom_line(aes(y = cumsum(rDHI_acum)/k), group = 1, col="blue") +
    scale_y_continuous(sec.axis = sec_axis(trans = ~ . / (1/k),
                                           name = " LLuvia acumulada [mm]", 
                                           breaks = seq(min(cumsum(x$rDHI_acum)),
                                                        max(cumsum(x$rDHI_acum)),
                                                        by=50)),
                       breaks = seq(min(x$rDHI_acum),
                                    max(x$rDHI_acum),
                                    by=10))+
    ggtitle(titulo) +
    theme(plot.title = element_text(hjust = 0.5))+
    geom_line(aes(y=aport_mean_SMA_diario/k), col="green",alpha= 0.9)

}, titulo=names(Daiyli_avgs4))

## $`-8.04904174804688__42.5383224487305`

## 
## $`-7.31787109375__42.5615844726562`

## 
## $`-8.05770874023438__42.6729469299316`

6.1 comprobación correlacion diferencia de nivel y aportacion

Realizamos una gráfica de la diferencia de nivel respecto a la aportación.

x<- cut_train[[1]]
plot(diff(x$nivel_mean_SMA)*20,  
     ylim = c(-1.5,2), 
     type = "l", 
     xlab = "",
     xaxt="n",
     yaxt="n",
     ylab = " ")
lines(x$aport_mean_SMA/100-1, col="red")

Vemos que existe un disparo en la gráfica… eso es porque existe un salto temporal de 3 días en diciembre… en los que el nivel vario 80 cm… cortamos los datos para eliminar ese disparo.

xx<- x[which(x$Date> ymd("2018/12/10")),]  
plot(diff(xx$nivel_mean_SMA)*20,  
     ylim = c(-1.5,2), 
     type = "l", 
     xlab = "",
     xaxt="n",
     yaxt="n",
     ylab = " ")
lines(xx$aport_mean_SMA/100-1, col="red")

Comprobamos las cross-correlation de estas dos variables.

ccf_nivel_aport<- ccf(diff(xx$nivel_mean_SMA),
                      xx$aport_mean_SMA,
                      lag.max = 1000)

print(paste0("Correlación máxima: ",
             max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
             ccf_nivel_aport$lag[
               which.max(ccf_nivel_aport$acf)], " horas"))

## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"

La aportación se ve reflejada en el nivel con un retraso de 15 horas con una correlación bastante alta de casi el 80%.

vec_cor<- vector()

for (i in seq(6, 144, by=2)) {
  x<- cut_train[[1]]
  x$prep_hourlySMA<- SMA(x$prep_hourly, i)
  x<- x[complete.cases(x),]

  rain_aport<- ccf(x$prep_hourlySMA, 
      x$aport_mean_SMA, 
      lag.max = 1000, 
      plot = F)
  
  vec_cor[i]<- max(rain_aport$acf)
  
  
}

plot(vec_cor, 
     xlab = "SMA points",
     ylab = "Corr")

Se ha decidido que el SMA ha emplear será de 48 puntos. A continuación realizamos una predicción de la aportación empleando diferentes métodos y variables.

• Aportación de hace 3 días
• Lluvia de hace 1 día
• Lluvia actual

path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r) r[which(r$Date > ymd("2018/10/01")),])
})


fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/12/01")
    
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})


cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
  y<- lapply(x, function(r){
    
    
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
}, 
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA")]
  
  
})


x<- cut_train[[1]]
x$aport1<- lag(x$aport_mean_SMA, 24)
x$aport2<- lag(x$aport_mean_SMA, 48)
x$aport3<- lag(x$aport_mean_SMA, 72)
x$aport4<- lag(x$aport_mean_SMA, 96)
x$aport5<- lag(x$aport_mean_SMA, 120)

x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)

x<- x[complete.cases(x),]


y<- cut_predict[[1]]
y$aport1<- lag(y$aport_mean_SMA, 24)
y$aport2<- lag(y$aport_mean_SMA, 48)
y$aport3<- lag(y$aport_mean_SMA, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)

y<- y[complete.cases(y),]

indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")

y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]




fit_1  <- svm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2  <- svm(aport_mean_SMA  ~ prep_hourlySMA48_lag + aport3 +  prep_hourlySMA48, data = x)
fit_3  <- svm(aport_mean_SMA  ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)

fit_4  <- lm(aport_mean_SMA ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5  <- lm(aport_mean_SMA  ~ prep_hourlySMA48_lag + aport3 + prep_hourlySMA48, data = x)
fit_6  <- lm(aport_mean_SMA  ~ prep_hourlySMA48_lag * aport3 * prep_hourlySMA48, data = x)


uncorrected<-y$aport_mean_SMA
pred_aport1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))

pred_aport4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_aport6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        aport3= y$aport3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))


plot(uncorrected, x=y$Date, type = "l",
     ylim = c(0,600), 
     xlab = paste0(range(y$Date)),
     ylab = "Aportación [m³/s]",
     main =  "Predicción de aportación \n empleando diferentes métodos")
lines(pred_aport1,x=y$Date, col="red")
lines(pred_aport2,x=y$Date, col="green")
lines(pred_aport3,x=y$Date,col="blue")
lines(pred_aport4,x=y$Date, col="red", lty=2)
lines(pred_aport5,x=y$Date, col="green", lty=2)
lines(pred_aport6,x=y$Date,col="blue", lty=2)

Del anterior gráfico podemos concluir que no vale la pena realizar la predicción aplicando svm (Support Vector Machine), o por lo menos no de esta forma…. Es interesante remarcar que solamente estamos en disposición de realizar un pronóstico a 3 días… dado que el modelo necesita el dato de aportación de hace 3 días para funcionar.

ccf(diff(x$nivel_mean_SMA), x$prep_hourlySMA48, lag.max = 1000)

print(paste0("Correlación máxima: ",
             max(ccf_nivel_aport$acf), " Se produce para un retraso de ",
             ccf_nivel_aport$lag[
               which.max(ccf_nivel_aport$acf)], " horas"))

## [1] "Correlación máxima: 0.79527611153213 Se produce para un retraso de -15 horas"

No tenemos valores de aportación en la página web, solamente nivel… por lo tanto hay que conseguir una buena correlación entre nivel y aportación. O directamente predecir el nivel con la lluvia.

path_hist_WRF<- here::here('Data/Parques/Belesar/Historico/Hist_D1D2_DHI_MERGED.RDS')
clean_data<- readRDS(path_hist_WRF)
clean_data_oct<-lapply(clean_data, function(x){
  y<- lapply(x, function(r){ 
    r<- r[which(r$Date > ymd("2018/10/01")),]
  r$nivel_mean_SMA1<- c(diff(r$nivel_mean_SMA),0)
  return(r)
  })
})


fecha_cut<- ymd("2019/01/15")
#cortar en entrenamiento y predicción
cut_train<- lapply(clean_data_oct[c(14,18,19)], function(x,fecha_end){
  y<- lapply(x, function(r){
    
    fecha_ini<- ymd("2018/12/01")
    
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
},
fecha_end=fecha_cut)
cut_train<- lapply(cut_train, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
  
  
})


cut_predict<- lapply(clean_data_oct[c(14,18,19)], function(x, fecha_ini){
  y<- lapply(x, function(r){
    
    
    fecha_end<- ymd("2019/02/20")
    Jan_data<- r[which(r$Date<fecha_end & r$Date>fecha_ini),]
    return(Jan_data)
    
  })
  return(y)
}, 
fecha_ini=fecha_cut)
cut_predict<- lapply(cut_predict, function(x){
  x[[1]][,c("Date","LON", "LAT", "prep_hourly", "Lluvia_mm", 
            "aport_mean_SMA", "nivel_mean_SMA","nivel_mean_SMA1")]
  
  
})


x<- cut_train[[1]]
x$nivel1<- lag(x$nivel_mean_SMA1, 24)
x$nivel2<- lag(x$nivel_mean_SMA1, 48)
x$nivel3<- lag(x$nivel_mean_SMA1, 72)
x$nivel4<- lag(x$nivel_mean_SMA1, 96)
x$nivel5<- lag(x$nivel_mean_SMA1, 120)

x$prep_hourly1<- lead(x$prep_hourly, 24)
x$prep_hourlySMA12<- SMA(x$prep_hourly, 12)
x$prep_hourlySMA24<- SMA(x$prep_hourly, 24)
x$prep_hourlySMA36<- SMA(x$prep_hourly, 36)
x$prep_hourlySMA48<- SMA(x$prep_hourly, 48)
x$prep_hourlySMA48_lag<- lag(x$prep_hourlySMA48, 24)

x<- x[complete.cases(x),]


y<- cut_predict[[1]]
y$nivel1<- lag(y$nivel_mean_SMA1, 24)
y$nivel2<- lag(y$nivel_mean_SMA1, 48)
y$nivel3<- lag(y$nivel_mean_SMA1, 72)
y$prep_hourly1<- lead(y$prep_hourly, 24)
y$prep_hourlySMA12<- SMA(y$prep_hourly, 12)
y$prep_hourlySMA24<- SMA(y$prep_hourly, 24)
y$prep_hourlySMA36<- SMA(y$prep_hourly, 36)
y$prep_hourlySMA48<- SMA(y$prep_hourly, 48)
y$prep_hourlySMA48_lag<- lag(y$prep_hourlySMA48, 24)

y<- y[complete.cases(y),]

indexxx<- 285
fecha_ini<- y$Date[indexxx]
fecha_end<- y$Date[indexxx]+ as.difftime(3, units = "days")

y<- y[which(y$Date> fecha_ini & y$Date< fecha_end), ]




fit_1  <- svm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48 , data = x)
fit_2  <- svm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag + nivel3 +  prep_hourlySMA48, data = x)
fit_3  <- svm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)

fit_4  <- lm(nivel_mean_SMA1 ~ prep_hourlySMA48_lag * prep_hourlySMA48, data = x)
fit_5  <- lm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag + nivel3 + prep_hourlySMA48, data = x)
fit_6  <- lm(nivel_mean_SMA1  ~ prep_hourlySMA48_lag * nivel3 * prep_hourlySMA48, data = x)


uncorrected<-y$nivel_mean_SMA1
pred_nivel1<- predict(fit_1, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel2<- predict(fit_2, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel3<- predict(fit_3, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))

pred_nivel4<- predict(fit_4, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel5<- predict(fit_5, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))
pred_nivel6<- predict(fit_6, data.frame(prep_hourlySMA48_lag =y$prep_hourlySMA48_lag,
                                        nivel3= y$nivel3,
                                        prep_hourlySMA48 =y$prep_hourlySMA48))


plot(uncorrected, x=y$Date, type = "l",ylim=c(0, 0.2),
     xlab = paste0(range(y$Date)),
     ylab = "nivelación [m³/s]",
     main =  "Predicción de nivelación \n empleando diferentes métodos")
lines(pred_nivel1,x=y$Date, col="red")
lines(pred_nivel2,x=y$Date, col="green")
lines(pred_nivel3,x=y$Date,col="blue")
lines(pred_nivel4,x=y$Date, col="red", lty=2)
lines(pred_nivel5,x=y$Date, col="green", lty=2)
lines(pred_nivel6,x=y$Date,col="blue", lty=2)